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Abstract. The level of current understanding of the physics of time-dependent strongly corre- 
lated quantum systems is far from complete, principally due to the lack of effective controlled ap- 
proaches. Recently, there has been progress in the development of approaches for one-dimensional 
systems. We describe recent developments in the construction of numerical schemes for general 
(one-dimensional) Hamiltonians: in particular, schemes based on exact diagonalization techniques 
and on the density matrix renormalization group method (DMRG). We present preliminary results 
for spinless fermions with nearest-neighbor-interaction and investigate their accuracy by comparing 
with exact results. 



1. INTRODUCTION 

The time evolution of strongly interacting quantum many body systems is one of the 
most challenging experimental and theoretical problems in physics. Recently, there has 
been progress in the investigation of the non-equilibrium time evolution of quantum 
systems in optical lattices [ 1] which has spurred theoretical interest in the time evolution 
of many body systems. The analytical treatment of such systems in equilibrium is 
restricted to only a few models, but in low dimensions efficient numerically exact 
methods have been developed and successfully applied to a variety of models in recent 
years. 

The perturbative Keldysh formalism can be applied to investigate the non-equilibrium 
behavior of such systems analytically [2]. Also, nonperturbative approximation schemes 
basing on functional integral techniques are under development in the context of high- 
energy physics yfl. Due to their perturbative or approximate character, these approaches 
cannot treat systems for arbitrary parameter values. Therefore, development of numeri- 
cal tools which do not have such limitations is crucial. In this contribution, we will focus 
on one-dimensional systems because the techniques developed in the past decade work 
very well for them, whereas it is rather difficult to apply these methods with the same 
accuracy and efficiency to higher dimensional systems. 

A number of different numerical methods for the investigation of one-dimensional 
strongly correlated quantum systems exist. The most important methods are Quantum 
Monte-Carlo (QMC) 0, Eh, exact diagonalization (ED) [6], and the density-matrix 
renormalization group (DMRG) L3.uil2l.lJI3]- Some approaches to time evolution using 
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QMC are known |Tl], [Tj] and others are under development B13I1 . but the numerical 



calculations are difficult to control. Therefore, we focus on the exact diagonalization and 
on the DMRG. Although it is possible to calculate all desired quantities using the ED 
for small system sizes, it is necessary to use the more involved DMRG method in order 
to reach sufficiently large system sizes to carry out a controlled finite-size scaling to the 
thermodynamic limit or to describe experiments with a large but finite number of sites, 
which are, e.g., realizable in optical lattices. In this contribution, we present results for 
the time evolution of a one-dimensional system obtained using an ED approach and use 
them for testing different DMRG approaches to the time evolution of the same system. 

The main difficulty in calculating the time evolution using the DMRG is that the 
effective basis determined at the beginning of the time evolution is not able, in general, 
to represent the state well at later times II 1411 because it covers a subspace of the system's 
total Hilbert space which is not appropriate to properly represent the state at the next 
time step. As will be shown in Sec. 13 it is possible that the representation of the time- 
dependent wave function very soon becomes quite bad. It is necessary either to mix 
all time steps into the density-matrix 1141 LUfl , or to adapt the density matrix. 

An approach for adaptive time-evolution basing on the Trotter- Suzuki decomposition of 
the time-evolution operator was developed in Refs. (ISO, EH- However, the method 
is restricted to systems with nearest-neighbor terms in the Hamiltonian only. In our 
approach, we use a general scheme closely related to exact diagonalization techniques to 
obtain the next time step during the time evolution, combined with an adaption scheme 
recently proposed by White ifloL l2(ill . Using this more general approach, it is possible 
to treat Hamiltonians containing arbitrary terms as long as the system can be efficiently 
handled with standard DMRG procedures. 

The article is organized as follows. First, we briefly review the DMRG method in Sec. 
El In Sec. we discuss general schemes for obtaining the time evolution of arbitrary 
quantum systems. In Sec.|4] we focus on strongly correlated systems and present results 
calculated using an exact diagonalization scheme. Subsequently, in Sec. 13 we discuss 
approaches and preliminary results obtained with non-adaptive and adaptive time evo- 
lution using the DMRG method and compare the results with the outcomes of the exact 
diagonalization approach. Finally, we discuss the results obtained and future work in 
Sec.|6l 

2. THE DMRG METHOD 

The DMRG method 0. Hi is described in detail in another contribution in this volume 
lElll . However, in order to better understand the main difficulties related to the calcu- 
lation of the time evolution of an initial state | y/o) using DMRG, we will give a short 
description of the main steps of the method here. The basic idea of the density-matrix 
renormalization group method is to represent one or more pure states of a finite sys- 
tem approximately by dividing the system in two and retaining only the m most highly 
weighted eigenstates of the reduced density matrix of the partial system. In combination 
with the numerical renormalization group approach (NRG) developed by Wilson [22] 
and the superblock algorithms developed by White and Noack Il23ll . this leads to a very 
powerful and efficient tool for the investigation of one-dimensional strongly correlated 
quantum systems on a lattice. 



RG-step 1 : Increase number of degrees of freedom: 
Add exact site to old system block 



Obtain |*F> 



oooo»«oooo 



Obtain p 
Diagonalize p 



New basis: eigenstates of p 
RG-step 2: Decrease number of degrees of freedom: 
Cutoff after m states 



Transform system block 
into new basis with only m states 



FIGURE 1. Sketch of the lattice and the flowchart of the DMRG iteration scheme. The left part of the 
lattice is the subsystem to which a site is added; the "sweep" proceeds from left to right. The flowchart 
shows the relevant steps of the DMRG procedure. For further details of the method see 1 21 ]. 



As depicted in Fig. I, the key steps are to increase the number of degrees of freedom 
of the partial system by adding sites, then to decrease the number of degrees freedom 
by retaining states below a cutoff. In this way, the method carries out a renormalization 
group procedure closely related to Wilson's NRG. 

In the first step of the algorithm, in which a site is added to one of the subsystems, the 
site-Hamiltonian is represented exactly in the site's many-bodybasis. This basic feature 
is exploited by the Trotter- Suzuki approach to time evolution 11171 flRfl to efficiently apply 
the local time-evolution operator in this part of the lattice. The subsystem's Hamiltonian 
is usually represented in an efficient reduced basis built up from the m most important 
eigenstates of the subsystem's reduced density-matrix. In the second step, the states 
one is interested in are obtained. These states are called "target states". In the original 
ground- state algorithm, these are the ground state and the few lowest lying excited states 
of the system, which are obtained by diagonalizing the Hamiltonian of the total system, 
e.g. by using the Lanczos diagonalization algorithm described in Sec. |4j However, one 
is not restricted to eigenstates of the Hamiltonian, and in this step may obtain alternative 
states. These properties are crucial for the time-evolution algorithms. The main problem 
is that the time-evolved state is not given by solving an eigenvalue problem. 

In the third step, the new effective basis is obtained by diagonalizig the reduced 
density matrix of the extended subsystem given by 



where the sum goes over all target states. In step four, only the m eigenstates with the 
largest eigenvalues are kept. The operators needed to represent the subsystem's Hamil- 
tonian, to form the pieces of the Hamiltonian connecting subsystems, and to calculate 
observables are transformed into this new reduced basis. This effective Hamiltonian of 
the subsystem is now the starting point for step one of the next iteration. In this way, 
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FIGURE 2. Comparison of the time evolution of the local density {«,) using full diagonalization and 
the Lanczos-approach described in Sec.0]for a system of spinless fermions on 14 lattice sites with open 
boundary conditions and interaction V (t = 0) = 0.5, V(t > 0) = 2.5. mi = 5 Lanczos vectors were kept. 
As can be seen, the error is negligible even after long times. 



every iteration step improves the accuracy of the obtained eigenstates and energies by 
improving the reduced basis used for the representation of the target states. 

Since both the Hamiltonian and the wavefunction | at ^ me t are represented in 
an incomplete basis, the result for the next time step \\j/(t + dt)} will have additional 
errors because the reduced basis is not an optimum representation for this state. In 
order to minimize these errors, it is necessary to form a density matrix whose m most 
important eigenvectors are "optimal" for the representation of the state \y(t)), as well 
as for \\ff(t + dt)) in the reduced Hilbert space. In Sec. 13 we will discuss ways of dealing 
with this problem. 



3. APPROACHES TO TIME EVOLUTION 

The dynamics of a quantum system is given by the solution of the time-dependent 
Schrodinger equation 

mj- t \ ¥ (t)) = H\ ¥ (t)) . (2) 

This equation is a first-order ordinary differential equation (ODE) which can be nu- 
merically solved directly by approximate integration schemes such as the Runge-Kutta 
method W24I1 . This method has been used within a non-adaptive DMRG scheme in Ref. 
ll25ll to obtain the time-dependent tunneling currents through a quantum dot connected 
to both non-interacting and interacting leads. This approach is a standard approach with 
an error e (dt) 4 , but does not conserve unitarity. The crucial numerical step in this pro- 
cedure is the H(t) I which must be implemented efficiently. An alternative implicit 
integration scheme conserving unitarity is the Crank-Nicholson procedure 12411 



(3) 



This approach has been used in Ref. within a non-adaptive DMRG approach to 
obtain the time evolution of a Bose-Hubbard system with an instantaneous change in 
the interaction strength at the beginning of the time evolution. In this method, the most 
costly numerical step is the calculation of the inverse (the denominator in Eq. O, which 
can be carried out using, e.g., a biconjugate gradient approach [18]. The accuracy of this 
operation also determines the error of this approach. Other, more involved implicit and 
explicit integration schemes are known (see, e.g., Ref. 92m ). but none of them has been 
used yet within the DMRG. 

An alternative to the direct integration of the Schrodinger equation is to treat the 
formal solution 

\ ¥ (t)) = U\ ¥ (t)} (4) 

directly, where U = e~ lHt / h is the time evolution operator. When full diagonalization of 
the Hamiltonian is viable, the time evolution can be expressed in the eigenbasis by 

W{t)) = Y,e- iEnt,h \n){n\w), 

n 

where E n and \n) are the eigenvalues and eigenfunctions of H and | y/o) is the initial state. 
However, for strongly interacting quantum systems it is not possible to fully diagonalize 
H for all but the smallest system sizes, so that an approximation to the time-evolution 
operator is needed. In the next section we present an ED approach with which such an 
approximation can be obtained within the Krylov basis to very high accuracy. 

4. TIME EVOLUTION USING EXACT DIAGONALIZATION 

Efficient iterative eigensolvers exist which can calculate the ground state and low ly- 
ing eigenstates of systems with Hamiltonians which can be represented as sufficiently 
sparse matrices. The full spectrum, however, cannot efficiently be obtained with these 
techniques. One example is the Ldnczos -procedure JzSEzIl which is presented in more 
detail in another contribution in this volume. In this procedure, the vectors of the Krylov 
subspace, the subspace spanned by vectors 

{\u },H\u },H 2 \u },...,H n \u }} , 

are orthogonalized with respect to the previous two vectors of the set, leading to the 
recursion relation 

\uj+i) = H\uj) - OCj\uj) - Pj\uj-i) (5) 

with the coefficients 



The Hamiltonian is then represented by the tridiagonal matrix 
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which can be easily diagonalized. For a review, see Refs. JES. If n is equal to the 
dimension of the Hilbert space, the approach is equivalent to a full diagonalization of H 
(albeit not a useful one due to numerical instability). However, it turns out that for the 
calculation of the ground state it is sufficient in most cases to carry out of the order of 
100 such iterations. 

It is possible to use a Krylov-space approach to approximate the time-evolution 
operator = e -' dt / hH from the time step t to the time step t + dt. In the following, 
we assume that the Hamiltonian H has no explicit time-dependence for t > 0. 

The time evolution through one interval IJ^t \ ty) can be approximated by J28ll29ll 



\w{t+dt)) = e- idt/hn \\j/(t)) w \ n {t) g-^W) \ T n (t) \\j/{t)) 



(8) 



where V n is the matrix containing all the Lanczos vectors \uj). The error in this approx- 
imation is given by (30] 
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Here || • || represents the euclidean norm and p = l^max — ^minl is the width of the 
spectrum of the Hamiltonian. This means that the number of Lanczos vectors and the 
size of the time interval dt needed to obtain a given accuracy depend on the total energy 
of the system. In general, the larger the system, the more Lanczos vectors are needed, 
assuming that dt and the filling are kept fixed. However, this is not a serious limitation 
because the characteristic oscillations take place on timescales of the order T c h ar ~ §• 
In order to fully resolve the dynamics it is necessary to adjust dt to be of the order of 
T c har anyways. Together with relation dTOl. for typical situations where an accuracy of, 
e.g., tol < 10 , is required, one finds that n < 20 is sufficient. Therefore, the matrices 
in Eq. ® are very small. The most time-consuming part of the calculation is then the 
multiplication H\u n ) needed to carry out the recursion to obtain the Lanczos vectors, it 
is important to implement this efficiently. Putting these steps together, the full procedure 
reads: 



1 . Estimate the number of Lanczos vectors needed for the given time step to obtain 
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FIGURE 3. Plot of the local density («,•) and the momentum distribution (n^) for a system of 18 sites 
with periodic boundary conditions, V(t = 0) = 0.5, V(t > 0) = 100 calculated with the ED approach 
presented in Sec. 0] One can clearly observe a collapse and revival of the discontinuity at the Fermi edge, 
while («,-) = for all times. In this calculation, dt = 0.001 and mi = 10 Lanczos vectors. 



2. Obtain V mL and T mL by performing the Lanczos iteration scheme with as 
starting vector. 

3. Compute \\]/(t + dt)} = Y n (t) e -idt/m„(t) y r^ \y(t)). 

4. Calculate observables. 

5. Continue starting with step 2 and replacing \y/(t)) by \\jf(t + dt)) until t mSLX is 
reached. 

Using this approach, we calculate the time evolution of a system of spinless fermions 
given by the Hamiltonian 

# = ~ l L ( c ]+ 1 c j + h ■ c - ) + v E n J n J+ 1 ' ( 1 1 } 

where we change the magnitude of the interaction strength V at t=0. The half-filled 
system is known to undergo a phase transition at the critical parameter value V = 2 
from a metallic (V < 2) to a CDW insulating phase (V > 2) 113 ill in the thermodynamic 
limit. Thus, by starting with a ground state obtained for V < 2 and applying the time- 
evolution operator with an interaction strength V > 2, one would expect oscillations in 
time between these two different phases. In addition, one would expect time-dependent 
Friedel-like oscillations in the local density (ni) near the boundaries for finite systems 
with open boundary conditions. Such oscillations are shown in Fig. |3 for two different 
times. In order to test the accuracy of the Lanczos calculations, we have compared to 
results obtained from full diagonalization of the Hamiltonian. As can be seen, the error 
remains negligible up to the fairly long times treated here. 

In Fig. results are shown for periodic boundary conditions for L = 18, a system 
size not accessible to full diagonalization. As can be seen, there is no change in the local 
density, for all times. This is expected due to translation symmetry, and is reproduced by 
the numerical calculations, which do not explicitly enforce the symmetry. 
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FIGURE 4. Flowcharts of adaptive time-evolution schemes. On the left, the procedure used in the 
Trotter- approach developed in lfl7t [T3l is sketched. On the right, the approach presented by White in 
I19l 12011 is shown. Note that the wave-function transformation is needed in both approaches (see Refs. 

mm)- 



At the initial time step at temperature T = 0, for V < 2, the momentum distribution 
(n(k)) has a singularity at the Fermi edge in the thermodynamic limit ll32ll . However, 
due to the limited system size, a step at kf is observed rather than a singularity. This 
discontinuity vanishes and then reappears in the course of the time evolution. This may 
be interpreted as a collapse and revival of the metallic state of the system. Note, however, 
that a step at the Fermi edge in a finite one-dimensional system does not automatically 
induce a metallic state in the thermodynamic limit. At present, it is unclear whether the 
singularity in the thermodynamic limit is revived in the course of the time evolution. 
This issue is under investigation and will be presented elsewherel33ll. 

Within full diagonalization, sizes L < 14 are possible for this system. For ED, L = 26 
can be reached on a supercomputer using a very basic code. By parallelizing the code and 
exploiting symmetries more completely, larger system sizes could be reached. However, 
they would nevertheless be small compared to the system sizes that we expect can be 
treated with the time-dependent DMRG. 



5. ADAPTIVE TIME EVOLUTION USING THE DMRG 

As pointed out in Sec.|2l the main problem in doing time evolution with DMRG is the 
update of the density matrix basis. In Fig.H we outline the algorithms for two different 
approaches to adaptive time evolution within the DMRG. In the first approach, a Trotter- 
Suzuki decomposition [34] of the time-evolution operator is used. The key feature of this 
method H 1 VL 1 1 8fl is the application of the exact local bond time-evolution operator to the 
two exactly treated sites for a particular superblock configuration, represented in Fig. [l] 
by filled circles. Results obtained by using this approach have been presented in Refs. 
lfl7l[l8Ll35ll3f3l . an extensive error analysis is given in J35ll . This approach is restricted 
to systems with nearest-neighbor terms in the Hamiltonian only because it relies on the 
fact that the sites to which the bond operator is applied are represented exactly. In order 
to formulate a more general method, it is useful to reexamine the approaches presented 
in Sec. |3j The most promising candidate is the Lanczos approach, presented in Sec. 
HJ which has small and well-controlled errors. An approach including states at all time- 
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FIGURE 5. Comparison between the time-evolution obtained with the ED approach of Sec.|4l DMRG 
without basis adaption, and DMRG with basis adaption, as described in Sec.|5]for a system of spinless 
fermions with 20 sites, open boundary conditions, and V(t = 0) = 0.5, V(t > 0) = 2.5, Ml = 10. 



steps in the density-matrix rather than basis-adaption has been used in Ref. [ 15] . We now 
describe a Lanczos-based approach using adaptive time evolution within the DMRG. 
The key ingredient is the implementation of the basis adaption in a way which represents 
the state \y(t)} as well as \ y(t + dt)) optimally. A scheme to do this has recently been 
proposed by White In this approach, the density-matrix basis is formed by 

including different time-steps within the time interval [t,t + dt]. As in standard DMRG, 
finite-system sweeps are performed until convergence is achieved. At this point, the 
density matrix basis has been optimized to represent both the state at time t and the state 
at the time t + dt. Observables at the time step t + dt can be calculated and | yf(t + dt) ) is 
then used as the starting point for propagation to the next time-step. In Fig.HJ we outline 
the procedure. While we currently include time steps t ,t + dt /3 ,t + 2dt/3, and t + dt 
in the density matrix, which intermediate states are optimal to include is still under 
investigation. At present, it is unclear how the sweeping influences the convergence 
behavior, although we believe that the accuracy of the state | \jf(t + dt)) can be improved 
by additional sweeping. Also, the errors in the reduced basis will accumulate with time. 
A more thorough error analysis and an investigation of possibilities to optimize this 
procedure are in preparation and will be presented elsewhere ll33ll . 

In Fig. we show preliminary results for a system of spinless fermions with open 
boundary conditions for L = 20. As can be seen, the time evolution without basis adap- 
tion yields a result that is qualitatively wrong even for very small times, while the adap- 
tive DMRG retains its accuracy for much longer times. A comparison of the non-basis- 
adapted and basis-adapted methods with our current time-evolution program, which does 
not utilize conserved quantum numbers W37I1 . shows that the error is significant even af- 
ter only 3 time steps for the non-adaptive method, while the adaptive method can reach 
50 time steps before the error becomes discernable. 



6. DISCUSSION 



In this contribution, we have proposed schemes for the time evolution of strongly 
correlated one-dimensional quantum systems based on exact diagonalization and on the 
DMRG. In the ED-based variant, we have applied a Krylov-space approach utilizing 
a Lanczos iterative diagonalization scheme to a system of spinless fermions. We find 
that an adiabatic change in the interaction strength induces a collapse and revival of 
the metallic state. The same system was used to explore possible methods for adaptive 
time evolution within the DMRG. As discussed in detail, it is crucial to adapt the 
density-matrix basis to the changing state at each time step; we have implemented such 
a scheme based on a proposal by White 119L l20ll . A comparison of results obtained 
with the adaptive DMRG scheme with the ED calculations shows that it can produce 
accurate results at fairly long time. However, a more thorough error analysis and a better 
understanding of how to optimize the scheme are necessary. Such work is in progress 
and will be presented elsewhere 13 311 . 
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